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ABSTRACT 

Hydrodynamical simulations of galaxy formation such as the Illustris simulations have pro¬ 
gressed to a state where they approximately reproduce the observed stellar mass function 
from high to low redshift. This in principle allows self-consistent models of reionization that 
exploit the accurate representation of the diffuse gas distribution together with the realistic 
growth of galaxies provided by these simulations, within a representative cosmological vol¬ 
ume. In this work, we apply and compare two radiative transfer algorithms implemented in 
a GPU-accelerated code to the 106.5 Mpc wide volume of Illustris in postprocessing in order 
to investigate the reionization transition predicted by this model. We find that the first gen¬ 
eration of galaxies formed by Illustris is just about able to reionize the universe by redshift 
z ~ 7, provided quite optimistic assumptions about the escape fraction and the resolution lim¬ 
itations are made. Our most optimistic model finds an optical depth of r ^ 0.065, which is 
in very good agreement with recent Planck 2015 determinations. Furthermore, we show that 
moment-based approaches for radiative transfer with the Ml closure give broadly consistent 
results with our angular-resolved radiative transfer scheme. In our favoured fiducial model, 
20% of the hydrogen is reionized by redshift z = 9.20, and this rapidly climbs to 80% by red¬ 
shift z = 6.92. It then takes until z = 6.24 before 99% of the hydrogen is ionized. On average, 
reionization proceeds ‘inside-out’ in our models, with a size distribution of reionized bubbles 
that progressively features regions of ever larger size while the abundance of small bubbles 
stays fairly constant. 

Key words: radiative transfer - methods: numerical - H II regions - galaxies: high-redshift 
- intergalactic medium - dark ages, reionization, first stars 


1 INTRODUCTION 

The cosmic microwave background (CMB) radiation was released 
when the Universe recombined at a redshift of z ~ 1100, leaving 
behind only a tiny residual free electron fraction. Yet in the present 
Universe, it is well established that the intergalactic medium (IGM) 
is highly ionized, as is inferred from the absence of Gunn & Pe¬ 
terson (1965) troughs in the absorption spectra of nearby quasars. 
Hence there must have been an ‘epoch of reionization’ (EoR) some¬ 
time in between, where photons emitted by stars and possibly 
quasars ionized the intergalactic hydrogen again (for reviews see 
Barkana & Loeb 2001; Fan et al. 2006a; Morales & Wyithe 2010). 
This is believed to first happen to hydrogen at z ~ 7 - 10, with he- 
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lium being ionized considerably later at z ~ 3 (Schaye et al. 2000; 
McQuinn et al. 2009) by the harder radiation of quasars. 

Observations of distant quasars show unambiguously that 
reionization was complete not much later than z = 6 (Fan et al. 
2006b). Observations of Lyman-cr emitters indicate rapid changes 
in their abundance at somewhat higher redshift (e.g. Ouchi et al. 
2010; Kashikawa et al. 2011; Caruana et al. 2014). Together with 
the relatively high optical depth to scattering on free electrons in¬ 
ferred from CMB observations (Bennett et al. 2013; Planck Collab¬ 
oration et al. 2014) and the discovery of early galaxies at z ~ 7 and 
higher (Bouwens et al. 2011; Pentericci et al. 2011; Oesch et al. 
2014), this suggests that reionization likely started considerably 
earlier than z ~ 7. The duration of the transition process, and the 
nature of the source population ultimately responsible for reioniza¬ 
tion, remain however subject of much observational and theoretical 
research. 
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Figure 1. Projected slices of thickness 7.1 cMpc through the binned gas and temperature fields of the 106.5 cMpc wide Illustris simulation at redshift z = 7, 
just before the externally applied UV background raises the temperature of the diffuse gas (the unit ‘cMpc’ stands for comoving Mpc). The left-hand panel 
shows the gas density field, with overlaid circles giving the locations of galaxies identified at this time by the group finder algorithm SUB FIND (Springel et al. 
2001). The right-hand panel gives the corresponding mass-weighted temperature field; the heated regions around the galaxies are caused by shocks associated 
with virialization and feedback-driven outflows. In turn, the high temperature there leads to collisional ionization of the gas. 


A particularly exciting prospect is that an observational break¬ 
through in this field may be imminent, in particular through a direct 
mapping of the EoR with 21-cm observations (see Zaroubi 2013, 
for a recent review). This has not yet been achieved, but impres¬ 
sive progress towards this goal has recently been made (e.g. Dillon 
et al. 2014; Parsons et al. 2014), and future instruments such as 
the Square Kilometer Array (SKA) or the James Webb Space Tele¬ 
scope (JWST) promise to revolutionize our understanding of the 
early universe, and of reionization in partieular. 

Numerous theoretical models for cosmic reionization have 
been constructed, often based on semi-analytic models of the reion¬ 
ization process or simple radiative transfer postprocessing of dark 
matter simulation outputs. Furlanetto et al. (2004) developed an ex¬ 
cursion set approach to reionization that has seen widespread use in 
analytic and semi-numerical models of reionization (e.g. Alvarez & 
Abel 2007; Zahn et al. 2011; Mesinger et al. 2011; Raicevic et al. 
2011; Battaglia et al. 2013). Many different numerical algorithms 
for direct radiative transfer simulations have also been developed 
over the years, in most cases however based on static density fields 
derived from dark matter only simulations or from simplified hy¬ 
drodynamic simulations (e.g. Sokasian et al. 2001; Ciardi et al. 
2003; Iliev et al. 2006; McQuinn et al. 2007; Zahn et al. 2007; 
Croft & Altay 2008; Trac et al. 2008; Aubert & Teyssier 2010; Ahn 
et al. 2012). 

Recently, full radiation hydrodynamics simulations that fol¬ 
low cosmic reionization and structure formation simultaneously 
and self-consistently have become possible. These calculations can 
in principle account for radiative feedback processes on forming 
galaxies, for example from inhomogeneous photoionization heat¬ 
ing. Pioneering work of this type has been presented by Gnedin 
(2000), but only in recent years it has become possible to study ap¬ 
proximately representative cosmological volumes in this way (e.g. 
Petkova & Springel 2011a; Paardekooper et al. 2013). 

Some of the most advanced studies of this kind include the 


simulations recently presented by Norman et al. (2015) and So et al. 
(2014), who use full radiation hydrodynamics simulations of cos¬ 
mic structure formation on a uniform grid. A fixed spatial grid reso¬ 
lution does however not allow a proper resolution of internal galaxy 
structure, which compromises the ability of the simulations to reli¬ 
ably predict the build up of stellar mass. To remedy this problem, 
Gnedin (2014) and Gnedin & Kaurov (2014) employ adaptive mesh 
refinement techniques and simulate galaxy formation in cosmologi¬ 
cal volumes with high spatial and mass resolution. Similar work has 
recently been presented by Pawlik et al. (2015), based on hydrody- 
namical SPH simulations coupled self-consistently to the radiative 
transfer scheme TRAPHIC (Pawlik & Schaye 2008). However nei¬ 
ther group evolves the simulations significantly past the EoR (let 
alone to z = 0) due to the high computational cost involved, so it 
is not yet clear whether these simulation models would also yield a 
plausible present-day galaxy population. 

This body of theoretical works has made it clear that the 
source population primarily responsible for reionization is most 
likely star-forming small galaxies at high redshift. While so-called 
Pop-III stars may boost the high redshift photon production rate, 
the overall contribution of these “first star” sources is likely only 
of secondary importance (Paardekooper et al. 2013; Wise et al. 
2014). Similarly, another potential source of ionizing photons, ac¬ 
tive galactic nuclei (AGN), is not expected to be critical at high 
redshift due to their large mean separation (Hopkins et al. 2007; 
Faucher-Giguere et al. 2009) and still fairly limited cumulative lu¬ 
minosity. Instead, it is often argued that small proto-galaxies, with 
stellar masses just around 10"^ M© or even lower may dominate 
the ionizing budget (Paardekooper et al. 2013). Ahn et al. (2012) 
showed that these mini halos alone can not complete reionization, 
but significantly contribute towards an earlier onset of reionization 
and thus enhance the optical depth towards the last scattering sur¬ 
face. These small halo mass systems are further assisted by sugges¬ 
tions that the escape fraction may strongly rise towards small halo 
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masses, as inferred by Wise et al. (2014) based on radiation hy¬ 
drodynamics simulation of faint high redshift galaxies. Using this 
finding in a semi-analytic model for reionization, Wise et al. (2014) 
have also demonstrated that the first galaxies may plausibly con¬ 
stitute the reionizing sources, yielding an optical depth consistent 
with Planck Collaboration et al. (2014) without exceeding the UV 
emissivity constraints by the Ly-cr forest. 

It is interesting however to note that the cosmic star formation 
history inferred from observations is predicted to rapidly decline 
towards high redshift (e.g. Ellis et al. 2013). This is also expected 
(Hernquist & Springel 2003) and desirable on theoretical grounds 
(Scannapieco et al. 2012), because simulation models of galaxy for¬ 
mation need to resort to an extremely efficient suppression of high- 
redshift star formation in order to successfully describe present-day 
galaxy properties (e.g. Stinson et al. 2013). A very low level of high 
redshift star formation is however quite the opposite of what seems 
necessary to explain early reionization and the comparatively high 
optical depth inferred from CMB experiments. This tension makes 
it difficult to attribute reionization entirely to young galaxies with 
more or less ordinary stellar populations. In ACDM models with 
low normalization erg, or in alternative warm dark matter cosmolo¬ 
gies, this problem is further exacerbated (Yoshida et al. 2003a,b), 
whereas in certain non-standard dark energy models that shift struc¬ 
ture growth to earlier times, for example ‘early dark energy’ (Wet- 
terich 2004; Gross! & Springel 2009), it may be alleviated. 

Recently, cosmological hydrodynamic simulations of galaxy 
formation such as the Illustris (Vogelsberger et al. 2014b) or Eagle 
(Schaye et al. 2015) projects have advanced to a state where they 
produce realistic galaxy populations simultaneously at z = 0 and at 
high redshifts, throughout a representative cosmological volume. 
This is achieved with coarse sub-resolution treatments of energetic 
feedback processes that regulate star formation, preventing it from 
exceeding the required low overall efficiency. Ideally the process of 
reionization should be coupled dynamically to such a galaxy for¬ 
mation model. However, running several such reionization simula¬ 
tions for different escape fraction parameterizations at the resolu¬ 
tion and size of Illustris would be computationally very expensive. 
Thus we have here reverted to study reionization in post process¬ 
ing only. This approach is not as self-consistent as the most recent 
generation of full radiation-hydrodynamics simulations of galaxy 
formation (Gnedin 2014; Norman et al. 2015; Pawlik et al. 2015), 
but it can take advantage of a successful description of the evolving 
source population and the intervening gas distribution down to low 
redshift. 

An important manifestation of feedback are galactic winds 
and outflows that substantially modify the distribution of the dif¬ 
fuse gas in the circumgalactic medium (CGM) and the IGM. This 
in turn also influences the gas clumping and the recombination in 
models of cosmic reionization. It thus becomes particularly inter¬ 
esting to test whether detailed models of galaxy growth such as 
Illustris are in principle also capable of delivering a successful de¬ 
scription of cosmic reionization, and if so, what assumptions are 
required to achieve such a success. 

This is exactly the goal of this paper. We use a sequence of 
snapshots with high time resolution of the high-resolution Illustris 
simulation and combine them with a radiative transfer scheme that 
is capable of accurately evolving ionizing radiation for an arbitrary 
number of sources. We are particularly interested in the question of 
whether the star formation history predicted by Illustris can reion¬ 
ize the universe early enough to be consistent with observational 
constraints, and how the reionization transition proceeds in detail 
in this scenario. Because we have implemented two different ra¬ 


diative transfer methods, we can also evaluate how well they inter¬ 
compare, thereby providing an estimate for systematic uncertain¬ 
ties related to these radiative transfer methods. We also explicitly 
test the impact and accuracy of the often adopted reduced speed of 
light approximation. 

This work is structured as follows. In Section 2, we introduce 
the Illustris simulation and our different radiative transfer schemes 
that we use to follow cosmic reionization. We then turn to a discus¬ 
sion of our primary results in Section 3. In Section 4, we analyse 
numerical caveats such as resolution dependence or systematic ef¬ 
fects due to different radiative transfer approximations and discuss 
our findings. This is followed by our conclusions in Section 5. 

2 METHODS 

2.1 The Illustris simulation 

Recently, Vogelsberger et al. (2014b) introduced the Illustris sim¬ 
ulation suite, an ambitious attempt to follow cosmological hydro¬ 
dynamics and the feedback processes associated with galaxy for¬ 
mation in a sizable region of the universe. The highest resolution 
simulation of the project employed 2 x 1820^ particles and cells in 
a box 106.5 Mpe across, yielding a mass resolution of 1.26x 10^ M© 
in the baryons, and 6.26 x 10^ M© in the dark matter. The cosmol¬ 
ogy adopted is given by = 0.2726, = 0.7274, = 0.0456, 

cTg = 0.809, Hs = 0.963, and Hq - 70.4kms"^Mpc"^ which is 
consistent with the most recent determinations from WMAP9 and 
Planck. The simulations employed the moving-mesh code AREPO 
(Springel 2010), which is well-suited to applications in cosmic 
structure formation. 

The physics model employed by Illustris includes radiative 
cooling, metal enrichment based on 9 elements, star formation, stel¬ 
lar evolution and mass return, supernova feedback by means of a 
kinetic wind feedback, and black hole growth and associated feed¬ 
back processes in a quasar- and radio-mode. We refer to Vogels¬ 
berger et al. (2013) and Torrey et al. (2014) for a description of the 
full details and basic tests of the model. A number of different stud¬ 
ies have analysed structure formation in Illustris, making it clear 
that many basic properties of the observed galaxy populations are 
approximately reproduced by the simulation model. This in partic¬ 
ular includes constraints on the stellar mass function at different 
epochs (Vogelsberger et al. 2014a; Genel et al. 2014), the mor¬ 
phologies and spectra of galaxies (Torrey et al. 2015), the colors 
of satellite systems (Sales et al. 2015), the stellar halos of galaxies 
(Pillepich et al. 2014), the nature of high redshift, compact galaxies 
(Wellons et al. 2015), the galaxy-galaxy merger rate (Rodriguez- 
Gomez et al. 2015), the kinematics and metal abundance of damped 
Lyman-alpha absorbers (Bird et al. 2015), or the evolution of the 
black hole mass density and quasar luminosity function (Sijacki 
et al. 2015). The galaxy formation predictions by Illustris are hence 
in broad agreement with observations, which adds additional moti¬ 
vation to ask whether they at the same time yield a plausible reion¬ 
ization history. 

A self-consistent treatment of the UV background using ra¬ 
diative transfer was however not included in Illustris, as this is 
still beyond reach in such large cosmological simulations that are 
evolved to low redshift. Instead, an external, spatially uniform and 
time-dependent UV background was imposed based on the model 
of Faucher-Giguere et al. (2009), and dense gas is self-shielded 
from UV background radiation using a prescription derived from 
Rahmati et al. (2013). A coarse treatment of an AGN proximity 
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effect was included where accreting AGNs modify the ionization 
balance of gas in their environment (Vogelsberger et al. 2013). In 
this work, we therefore aim to study reionization through postpro¬ 
cessing of the Illustris simulation, making use of the significant 
number of output dumps (more than 128, with an output spacing 
as in the Aquarius project, Springel et al. 2008) stored for the cal¬ 
culation. This allows us to take the temporal information about the 
growth of cosmic structures into account, avoiding the simplifica¬ 
tion of a static density field often adopted in past work. Another ad¬ 
vantage of Illustris is the reasonably large box size of 106.5 Mpc. 
While a still larger volume would clearly be desirable, studies of 
cosmic variance of reionization have suggested that ~ 100Mpc 
corresponds to the minimum box size required to obtain a reliable 
mean reionization history (Mesinger & Furlanetto 2007 ; Iliev et al. 
2014). The volume available in Illustris falls slightly short here, but 
is approximately still sufficient for hydrogen reionization. Studying 
Hell reionization with Illustris would be more problematic however 
due to the incomplete sampling of bright and rare quasars. 


2.2 Radiative transfer 


The governing equation of radiation transfer is the Boltzmann 
equation for the photon distribution function. 


dfy 

dt 


dx 




dfy 

dt 


dfy 

dt 


( 1 ) 


with cosmological scale factor a and speed of light c. The function 
fy{x, n, y, t) describes the number density of photons of frequency v 
at time t, propagating in the direction of the normal vector n. In the 
following, we will replace the frequency dependence by one effec¬ 
tive frequency bin that accounts for all hydrogen ionizing photons, 
for simplicity. This can in principle be straightforwardly general¬ 
ized to multiple frequency bins, as is necessary, for example, to also 
treat helium reionization. In equation (1), we ignore the redshifting 
of ionizing photons which is justified if photons are destroyed in 
ionization events shortly after their creation. 

Solving equation (1) in full generality is very difficult, hence 
a large variety of approximation methods for cosmological radia¬ 
tive transfer have been developed, each with different advantages 
and shortcomings that ultimately dictate the regimes where they 
are applicable. Many radiative transfer schemes are based on the 
idea of characteristics (Mihalas & Mihalas 1984), where the op¬ 
tical depth is individually integrated along rays between different 
computational cells. As this operation quickly becomes extremely 
costly many different techniques have been developed to accel¬ 
erate the calculations, for example by using short-characteristics 
(e.g. Nakamoto et al. 2001), or by hierarchically splitting rays 
(Abel & Wandelt 2002; Trac & Cen 2007). Ray-tracing schemes 
are particularly well suited for dealing with a small number of iso¬ 
lated sources, and in that sense are ideal for following the growth 
of an isolated ionized region around a point source (Abel et al. 
1999). However, the reionization problem involves a large number 
of sources (all the stars in a large number of galaxies), favouring 
the use of other methods. 

One popular approach is to take moments of the radiative 
transfer equation, leading to an evolution equation for the mean 
specific intensity that needs to be closed with an estimate of the 
local Eddington tensor (Gnedin & Abel 2001; Aubert & Teyssier 
2008; Petkova & Springel 2009; Finlator et al. 2009). The simplest 
variant of this approach is fiux-limited diffusion (ELD) of radia¬ 
tion (Levermore & Pomraning 1981; Turner & Stone 2001), which 
works particularly well in the optically thick regime but may fail 


badly in situations where the optical depth is low and shadowing 
might be important. Better accuracy can be obtained in moment- 
based approaches when the local Eddington tensor can be estimated 
with reasonable accuracy. This is attempted, for example, in the op¬ 
tically thin variable Eddington tensor (OTVET) approximation of 
Gnedin & Abel (2001), or by invoking other simple approxima¬ 
tions such as the so-called Ml-closure (Levermore 1984; Aubert & 
Teyssier 2008). 

More general radiative transfer such as Monte Carlo methods 
(e.g. Maselli et al. 2003; Nayakshin et al. 2009) may also be in¬ 
voked, but their computational cost is comparatively large making 
it difficult to avoid a high level of noise in large-scale simulated ra¬ 
diation fields. The TRAPHIC approach of Pawlik & Schaye (2008) 
improves on this by restricting the transport of photon packets to 
a finite set of angular cones, while retaining the flexibility of the 
Monte Carlo scheme. A conceptually similar idea is followed in the 
radiation advection scheme of Petkova & Springel (2011b), where 
the radiation field is directly discretized in angular space and trans¬ 
ported with a conservative advection solver on an unstructured grid. 

In our work here, we consider two different numerical meth¬ 
ods to solve the advection part of the radiation transport, described 
by the left hand side of equation (1). The first method is the cone- 
based advection method proposed by Petkova & Springel (2011b), 
and the second method is a simpler moment-based approach where 
the Ml closure for the Eddington tensor is used (Aubert & Teyssier 
2008; Rosdahl et al. 2013). This allows us to assess how strongly 
fundamental differences in the numerical radiative transfer approxi¬ 
mation affect the predicted reionization transition, thereby yielding 
an estimate for this contribution to the systematic uncertainties in 
theoretical EoR predictions. 


2.3 Cone-based advection method 


In the cone based method introduced by Petkova & Springel 
(2011b), the angular dependence of the distribution function fy is 
discretized in cones of equal opening angle. We use the HEALPIX 
tessellation scheme (Gorski et al. 2005) to subdivide the unit 
sphere, giving Apix = 12 X 4^ cones of equal solid angle, where 
n is an integer parameter determining the resolution level. The spa¬ 
tial dependence can be discretized using either a structured or un¬ 
structured mesh. As we are mostly interested in volume-weighted 
effects of reionization (such as the volume filling factor of ionized 
gas) and want to efficiently exploit GPUs as computational accel¬ 
erators we adopt a Cartesian grid in this work. This leaves us with 
Apix photon fields, // at grid coordinates Xq^^^ with fy -Y^fi 

being the total angle-integrated intensity. 

Leaving the source terms aside for the moment (which are 
treated in an operator-split fashion), the cone-based method solves 
a hyperbolic advection equation. 


dt dx 



= 0 , 


( 2 ) 


for each of the angular-decomposed photon fields f. The principal 
transport direction rii of cone / is set by the HEALPIX tessellation 
of the unit sphere. However, to make sure that the entire cone of 
each tile is illuminated homogeneously and that the transport does 
not excessively collimate the radiation into an angle smaller than 
the prescribed angular resolution, rii is obtained by taking 



( 3 ) 


and additionally, if this direction should fall outside of a cone with 
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center hi and opening angle 0max = >/4/A/pix, it is projected back 
inside that cone. Here hi is the geometric center of the HEALPIX 
cell. In other words, the advection direction is constrained to lie 
within the cone but otherwise tries to reduce gradients in the pho¬ 
ton intensity across the cone, thereby ensuring a homogenous illu¬ 
mination of the cone. 

We typically integrate the above transport equations using a 
reconstructed piecewise constant photon intensity field combined 
with an upwind scheme to determine the photon flux between cells. 
We use explicit time integration and hence need to restrict the 
timestep to something of order h/c, where h is the cell radius and 
c the speed of light. This leads to quite small timesteps. This could 
be mitigated by invoking a reduced speed of light approximation 
(Gnedin & Abel 2001), but as we discuss in section 4.3 we have 
not used this here to avoid any inaccuracies that can be introduced 
by this approximation. 


2.4 Moment method with Ml closure 


For comparison with the angular resolved transport scheme, we 
also implemented a moment based advection scheme. Here the ra¬ 
diative transfer equation is discretized by taking the first two mo¬ 
ments of equation (1): 


- + VF.0. 

^ + c^VP = 0, 
dt 


(4) 

(5) 


where the photon density N, the fiux vector F, and the photon 
pressure tensor P have been introduced. Solving this system of 
equations requires a closure for the radiation pressure tensor P. 
The simplest closure is to assume an isotropic pressure, which 
yields the family of fiux-limited diffusion methods (Levermore & 
Pomraning 1981). Another possibility is the OTVET approxima¬ 
tion where a preferred local streaming direction of photons is esti¬ 
mated based on the source locations and their strengths (Gnedin & 
Abel 2001). More accurate estimates of a variable Eddington ten¬ 
sor (VET) may for example be obtained by employing a special 
short-characteristics method to estimate the local Eddington tensor 
(Davis et al. 2012; Jiang et al. 2014). 

An alternative is the Ml method, where a fully local clo¬ 
sure is used instead (Levermore 1984; Aubert & Teyssier 2008). 
This greatly simplifies the computations because information about 
the surrounding regions does not explicitly have to be taken into 
account. The photon pressure tensor of the Ml closure can be 
parametrized as: 


P 

with 




-e 3^-1 

1+ — n^n\N, 


"■W\- 


3 + 4/" ^ 

5 + 2 V4 - 3/2 ’ cN' 


( 6 ) 


(7) 


and identity matrix I. Here / essentially determines how directed 
the fiux is, interpolating between the two limiting cases of radiation 
diffusion and photon streaming. The parameter ^ varies smoothly 
from I to 1 between these cases. In the case of ^ which oc¬ 
curs for an undirected flux, the result is isotropic advection. The 
other limit of ^ = 1 corresponds to fully directed transport with one 
dominant propagation direction. 

While in practice the Ml closure often produces surprisingly 
accurate results, it is easy to construct situations where it fails. Eor 


example, one obvious shortcoming occurs when two light beams 
directly encounter each other from opposing directions. As the pho¬ 
ton field is essentially described by Ml as a collisional fluid with 
a restricted form for the local Eddington tensor, a ‘scattering’ at 
the beam intersection point is unavoidable, resulting in an unphysi¬ 
cal propagation direction. Nevertheless, one can hope that in many 
practical applications such inaccuracies occur sufficiently rarely 
that the overall results are still reliable. Whether or not this is really 
the case is problem dependent and needs ultimately be tested by 
comparing with more accurate techniques that are based on differ¬ 
ent approximations, something that we also pursue in this work. 


2.5 Ionization network and time integration 


The evolution of the ionization state of hydrogen is given by the 
following system of equations: 

driy T . 

— = -ccrniiXiiy + [otaCT’) - aBCT")] «h (1 “ •*) . (8) 

- (^A{T)n]^ (1 - x)^ -l5{T)nl^x{\ - x) - ccrnuxny, (9) 

with the hydrogen number density , the frequency averaged ion¬ 
ization cross section cr, the number density of ionizing photons Uy 
and the neutral hydrogen fraction v. The case A and case B re¬ 
combination rates and and collisional ionization rate p are 
taken from Hui & Gnedin (1997). The chemical network is comple¬ 
mented by source terms in the thermal energy evolution describing 
cooling and photoionization heating: 


du 


= EyCO-riYiXny - C(r), 


( 10 ) 


where Sy is the average gain in thermal energy per ionization event. 
The cooling rate C{T) includes cooling due to recombination, col¬ 
lisional ionization, collisional excitation and Compton cooling by 
scattering on CMB photons. 

We solve the chemical network following the method outlined 
in Rosdahl et al. (2013). Eirst, the photon number Uy is updated 
by an implicit Euler step, then the temperature u and finally the 
ionization state v. Each partial update step is implicit in the already 
updated quantities and the currently updated one. The full radiative 
transfer equation (1) is then solved by Strang-like operator splitting, 
computing first the left hand side for half a time step (advection), 
then applying the source and sink terms, and finally completing the 
timestep by another half advection step. The two advection half 
steps themselves are solved by dimensional splitting. 

In our cone-based advection method, we usually use the rz = 1 
HEALPIX tessellation level with - 48 pixels. Due to the large 
number of photon fields the computations are demanding in mem¬ 
ory as well as in computational power. To speed up our simulations, 
we are employing GPUs, an approach similar in spirit to Aubert & 
Teyssier (2010). We use a uniform Cartesian mesh which makes 
our algorithm ideal for GPU computing, as this leads to a regular 
memory access pattern and facilitates the use of groups of threads 
with the same execution path. Details of our GPU implementation 
will be described elsewhere (Bauer et al., in preparation). 


2.6 Reionization in post processing 

In this work, we study the progress of cosmic reionization by fol¬ 
lowing the time-resolved evolution of the density field and the 
source population of the Illustris simulation in postprocessing. The 
unstructured Voronoi mesh that stores the baryon distribution in 
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the simulation outputs needs to be rebinned onto a regular Carte¬ 
sian mesh to allow use of our GPU-based code. To this end, we 
assign the mass of each Voronoi cell onto the density grid using 
a spline-kernel assignment. This yields a less noisy and smoother 
density field than obtained, e.g., by using a clouds-in-cells (CIC) or 
nearest-grid-point (NGP) assignment kernel (Hockney & Eastwood 
1981). The actual density field used in our reionization calculation 
is then continuously updated by linearly interpolating in time be¬ 
tween the two nearest binned density grids available in log a space, 
where a - 1/(1 -h z) is the cosmological scale factor. The Illustris 
outputs are spaced roughly 65 Myrs apart at the relevant redshift 
z ~ 6, allowing us to bin the density field in total 35 times between 
redshift z = 21.8 and z = 4.9. Using a set of small subboxes cut 
out from Illustris and stored with much higher time-resolution we 
have checked that the sparser time resolution of the main outputs is 
still reasonably accurate. It shifts the completion of reionization to¬ 
wards slightly earlier times, but the uncertainty in our results is still 
dominated by the parametrization of the escape fraction. Initially 
we start with a uniform temperature field with T = 100 K and then 
follow the temperature evolution using equation (10), ignoring the 
intrinsic temperature evolution of the underlying Illustris simula¬ 
tion. We note that in this paper we consider hydrogen reionization 
only. 

2.7 Escape fraction 

Of all ionizing photons emitted by stars, only a fraction /esc reaches 
the IGM, while the other photons are absorbed by the denser inter¬ 
stellar medium (ISM) or by dust. In principal, the escape fraction 
depends on individual halo properties like halo mass or dust con¬ 
tent. Unfortunately, only little is known about the real values of the 
escape fractions, especially at high redshifts, making this parameter 
one of the primary uncertainties in studying the EoR. However, the¬ 
oretical models have started to constrain the escape fraction, albeit 
with large systematic uncertainties. Eor example, the simulations 
of Wise et al. (2014) suggest a sharp increase of the escape fraction 
towards smaller masses, reaching 50% at halo masses of 10^ M©. 

The simplest model is evidently to adopt a globally constant 
escape fraction that is the same for all galaxies at a given epoch. 
This is what we shall assume here, because one may well argue 
that in light of the many other uncertainties adopting a more com¬ 
plicated model would not be justified (see also the discussion in 
Pawlik et al. 2015). As a default choice for the constant escape 
fraction model we have considered /esc = 0.2 (C2 model). This al¬ 
most certainly constitutes an overestimate for low redshift galaxies, 
but is perhaps not overly optimistic at high redshift. 

Besides such a globally constant escape fraction, we also con¬ 
sider an evolution of the /esc value with time. To this end, we use 
the model of Kuhlen & Eaucher-Giguere (2012) who proposed an 
evolution of the escape fraction as a function of redshift according 
to: 


Our default choices for the parameters /o and k are /o = 0.04 and 
K - A, implying a rise of the escape fraction with redshift (VI 
model). We have also calculated results for a variety of other fidu¬ 
cial parameter choices as well. Eor an overview see Table 1. 

We stress again that the escape fraction is essentially a free 
parameter in our treatment. Its interpretation is complicated by the 
fact that the absorbing ISM is not totally absent in our reionization 
simulations. It is just severely underesolved, and thus some part 



Figure 2. Time evolution of the total ionizing luminosity resulting for our 
different escape fraction models. The model with a constant escape frac¬ 
tion (blue lines) tracks the shape of the cosmic star formation history. Our 
default models with a redshift dependent escape fraction are given by the 
green and red lines. For comparison, we also include the ionizing luminos¬ 
ity implied by the quasars included in our simulation (dashed black line), 
adopting a simple conversion of quasar accretion rate to ionization radiation 
output. This demonstrates that the quasar contribution picks up too late to 
cause the initial hydrogen reionization, but it may play a role in keeping the 
universe ionized at later times, as well as for completing late-time helium 
reionization. 

of the boost in recombination rates due to a highly clumped envi¬ 
ronment is missing. However, photons are still consumed for (re¬ 
peatedly) ionizing the material in these high density regions. Unre¬ 
solved small galaxies and thus missing UV photons might be com¬ 
pensated by a slightly larger escape fraction than otherwise would 
be needed. 


2.8 Source modelling 

The main source of ionizing photons considered in our model is or¬ 
dinary stellar populations in young stars, which arguably appear to 
be the most likely source responsible for reionization. In this work 
we are mainly interested in testing this hypothesis based on directly 
adopting the stellar populations forming in the Illustris simulation. 
Their ionizing luminosities as a function of stellar age and metallic- 
ity are taken from STARBURST-99 (Leitherer et al. 1999). Eigure 1 
gives a visual impression of the binned gas density field in Illus¬ 
tris at z = 7, and of the clustered galaxy population that represents 
our source population. Because the local ionization rate can change 
on much shorter timescales than the density, we bin the luminosity 
field with much finer time resolution than the density field. In our 
default set-up, we binned it 200 times during the duration of the 
reionization simulation including the actual birth time of the stellar 
sources. The actual luminosity used in the simulation to integrate 
the source terms is then again interpolated in log a space from this 
large grid of luminosity density fields. The ionizing luminosity of 
a single stellar source is always distributed in a photon conserving 


MNRAS 000, 1-18 (2015) 










Hydrogen Reionization in Illustris 1 


Overview of radiative transfer simulation models 


Name Resolution /o k 

Advection method 


Vl_X_CONE 

256^ . 

..5123 

0.04 

4 

Cone based method 

V5_X_CONE 

256^ . 

..5123 

0.04 

3.6 

Cone based method 

C2_X_CONE 

256^ . 

..5123 

0.2 

0 

Cone based method 

V1_X_M1 

256^ .. 

. 10243 

0.04 

4 

Ml method 

V5_X_M1 

256^ .. 

. 10243 

0.04 

3.6 

Ml method 

C2_X_M1 

256^ .. 

. 10243 

0.2 

0 

Ml method 


Table 1. Overview of the different radiation transfer calculations performed for this study. We carried out runs with a resolution ranging from 256^ up to 
1024^ cells, which is indicated by replacing the placeholder W’ with the number of cells per dimension in the actual run name. We compare different escape 
fraction parameterizations, characterized by /o and k. The three different choices we adopted for this are labeled ‘VI’, ‘V5’ and ‘C2’ in the simulation names. 
Finally, for each of the models we compare two different advection schemes for the radiation, one based on the Ml closure relation, the other on an explicit 
discretization of the solid angle (‘cone based’). 


way to the radiation grids. The assignment of the luminosity is done 
using a CIC interpolation scheme. 

The ionizing sources we have at our disposal from the main 
Illustris simulation are related to the star formation rate and black 
hole accretion rate densities. The corresponding time evolution of 
the net ionizing luminosity density is shown in Figure 2, together 
with different scenarios for the escaping luminosity according to 
our escape fraction models. The blue line shows our constant es¬ 
cape fraction model (C2 model), while the green line is our default 
model for a time-varying escape fraction (VI model). The red line 
(V5 model) represents a variable escape fraction model with a dif¬ 
ferent value for the exponent k than in the V1 model. These models 
rapidly rise for some time, but at around a redshift of z - 8, the 
escaping radiation actually reaches a maximum and then starts to 
slowly decline again. Here the strong decrease of the escape frac¬ 
tion from 1 to 0.04 marginally over-compensates the further in¬ 
crease of the star formation density and the raw ionizing luminos¬ 
ity. We note that the maximum coincides with the epoch where we 
expect most of the hydrogen to be reionized. 

An alternative source of ionizing radiation is in principle pro¬ 
vided by AGNs. We assume a bolometric AGN luminosity de¬ 
scribed by 

= {l-€f)€rM^nc\ ( 12 ) 


with a radiative efficiency of - 0.2 and an energy fraction of 
(1 - 6f) = 0.95 available for radiation. The luminosity is converted 
into a rate of ionizing photons assuming a parameterized AGN SED 
(Korista et al. 1997) equal to 


yAGN(^) ^ ^auv 


exp 


hv 


exp 


10-^Ryd 

hv 


(13) 


with a suppression of the UV component at a temperature of 
7 "bb = 10^ K, a UV component slope of ctuv = “0-5 and an X-ray 
component slope of crx = -1. To obtain an approximation for the 
maximum contribution to reionization, we assume an escape frac¬ 
tion of /esc,AGN = 1- Comparison, we show the AGN ionizing 
luminosity as a dashed black line in Figure 2. The AGN contri¬ 
bution only becomes relevant at around z = 6, but by then most 
of the hydrogen must already be ionized, making it unlikely that 
AGNs are significantly contributing to the initial EoR transition. 
However, they might play a role in keeping the universe ionized 
at a later time, and almost certainly are important for completing 
helium reionization at lower redshift, thanks also to their harder 
spectrum. As we only study hydrogen reionization, in the following 
only stellar emission is considered as sources of ionizing photons. 


3 SIMULATION RESULTS 

An overview of the various reionization simulations carried out in 
this work is given in Table 1 . To examine resolution dependences, 
we computed several of our models with grid resolutions of up 
to 1024^ cells. We considered two different models with a time- 
variable escape fraction and contrasted them with one model with a 
fixed escape fraction. In order to asses the impact of different radia¬ 
tive transfer methods we performed most of our simulations both 
with our default cone-based approach as well as with the moment- 
based method with Ml closure. 

The progress of reionization in different environments is visu¬ 
ally shown in Figure 3, where we compare two regions around very 
massive halos with a more average environment around a typical 
medium-sized halo, and an underdense region. The different projec¬ 
tions show the four regions at six different output times. Reioniza¬ 
tion starts inside the most massive halos first, and then quickly ion¬ 
izes the surrounding regions. Compared to such a high-density en¬ 
vironment, the onset of reionization is considerably delayed around 
a medium mass halo. More drastically, the underdense region only 
begins to be reionized once the denser regions have almost com¬ 
pleted their reionization transition. The visual impression is thus 
qualitatively consistent with an inside-out reionization scenario in 
which halos in high-density regions are affected first and lower den¬ 
sity voids are reionized rather late, for the most part after overdense 
gas been reionized (Razoumov et al. 2002). This contrasts with sug¬ 
gestions that low density regions are reionized quite early and only 
then the reionization fronts progress to ionize filaments and gas in 
halos (Gnedin 2000). 


3.1 Reionization history 

In the upper panel of Figure 4, we show the volume weighted 
ionization fraction 1 - vy as a function of time for three of our 
high-resolution radiative transfer calculations. We see a rapid, ex¬ 
ponential rise at high redshift, and an approach to unity at around 
z - 6 - 8. The end of reionization and the rapid phase transition to 
a reionized universe is better visible in displaying the neutral frac¬ 
tion Ay, which is also included in the figure. For comparison, we 
also show observational constraints derived by Fan et al. (2006b) 
from quasar absorption lines (symbols with error bars). Especially 
our model V1_I024_M1 reproduces the suggested end of reioniza¬ 
tion in these data quite well, although the residual neutral fraction 
comes out slightly low. 

In the lower panel of Figure 4, we consider the evolution of 
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Figure 3. Progression of reionization as seen in the neutral hydrogen fraction in slices through selected sub-volumes in Illustris, each with a side-length of 
21.3 cMpc. Each column shows the time evolution of a different, randomly selected environment in our V1_1024_M1 model; the two columns on the left 
correspond to an average density higher than the mean, the other columns have medium and low mean density, as labeled. Each row gives a different redshift, 
from z = 9 (top) to z = 6 (bottom). 
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Figure 4. The ionization history for different escape fraction parameteri- 
zations. The upper panel shows the evolution of the volume weighted ion¬ 
ization fraction (solid lines rising to low redshift) and the neutral hydrogen 
fraction (lines dropping towards low redshift), as a function of redshift. For 
comparison, observational constraints by Fan et al. (2006b) are shown as 
symbols with error bars. The lower panel shows the ratio between mass- 
and volume-weighted ionization fractions (1 - xm) / (1 - ^v) for the differ¬ 
ent models. 

the ratio between mass- and volume-weighted ionization fractions, 
(1 - xm) / (1 - -^v)- This quantity is the average density of the ion¬ 
ized hydrogen compared to the average hydrogen density of the full 
box: 

1 — Xm ^ Vtot (1 “ -^^m) ^tot _ 1 ^ionized _ (p)ionized (14) 

1 “ -^V ^tot (1 “ -^v) l^tot (p)tot l^ionized (p)tot 

This ratio stays at or above unity for all time, which can be in¬ 
terpreted as a signature of an inside-out character of reionization 
(Iliev et al. 2006). Overdense environments around our sources ion¬ 
ize first, so that the ionized volume is always overdense on average. 
Interestingly, the evolution of the mean overdensity of ionized re¬ 
gions with time also differs for our different escape fraction models. 
The run assuming a constant escape fraction model starts reioniza¬ 
tion later but then progresses somewhat more rapidly. This model 
maintains the highest value of (1 - Xm) / (1 - -^v) for most of the 
simulated timespan. Here reionization is particularly biased to over- 
dense regions and is stuck there for a comparatively long time, until 
the final reionization transition occurs on a short timescale and the 
IGM at mean density is ionized as well. Interestingly, even though 
the variable escape fraction models show some variety in the time 
of the onset of reionization and the remaining neutral fraction, the 
evolution of (1 - Xm) / (1 - Xy) is still rather similar among these 
models. They begin reionization earlier and thus generally have low 
mean overdensities of the ionized volume at any given time. 

That the character of the reionization process is best described 


Figure 5. Time evolution of the neutral hydrogen fraction at given gas 
overdensities. Each of the lines shows the result for a small density range 
around a nominal comoving density threshold, as labeled, for our model 
V1_1024_M1. We see that the ionization of dense gas begins earlier, but 
this gas also ends up with a higher residual neutral fraction once reioniza¬ 
tion is completed. 

as an inside out transition can be seen in more detail in Figure 5, 
where the time evolution of the average neutral fraction is shown 
for regions of a fixed given overdensity. Highly overdense regions 
start to become ionized quite early on, assisted by collisional ion¬ 
ization in virialized halos. As a result, the reionization process for 
high density regions is also considerably less sudden than for lower 
density gas. After reionization is essentially completed, the behav¬ 
ior however reverses; now overdense regions show on average a 
final ionization degree that is considerably lower than for lower 
density regions. This can be understood as a result of the higher 
recombination rate in the denser regions, shifting the equilibrium 
value of the ionized fraction in a given UV background. Generally 
speaking, we find that denser regions start to ionize earlier but keep 
a higher neutral fraction than underdense regions. 

Reionization is clearly not an instantaneous transition but re¬ 
quires a certain amount of time. It is hence interesting to character¬ 
ize the epoch of reionization not only with a single redshift but also 
to ask how long the duration to a reionized universe takes. To this 
end we define the duration Az as the interval in redshift space dur¬ 
ing which the volume-weighted neutral fraction drops from 80% 
to 20%. Our fiducial model V1_1024_M1 leads to an extent of 
Az = 2.28 for the epoch of reionization, corresponding to a time- 
span of At = 251.2 Myr for this period. The V5_1024_M1 model 
shows a slightly longer duration of reionization with Az = 2.61 and 
At = 341.3 Myr. Reionization lasts only over a span of Az = 1.13 or 
At = 190.0 Myr in our C2_1024_M1 model. Thus our models with 
a variable escape fraction show a more extended epoch of reion¬ 
ization compared to the model with a constant escape fraction. In 
the models with a variable escape fraction, the ionizing luminos¬ 
ity is initially higher, but once most of the ionizing luminosity has 
become available, the variable escape fraction gradually begins to 
limit the amount of escaping UV radiation, resulting in a more pro¬ 
longed epoch of reionization. 

One can also ask whether galaxies of different stellar masses 
are all ionized at the same time, or whether there are significant 
systematic trends of the mean reionization epoch as a function of 
galaxy size. To this end, we have considered the sample of all z = 6 
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Figure 6. Mean reionization redshifts of the immediate surroundings of 
galaxies as a function of their z = 6 stellar masses for our highest reso¬ 
lution run of the VI and V5 models. The symbols give the means in each 
mass bin, with the error bars showing the statistical error of the mean. The 
solid lines illustrate the +lcr variance in each mass bin. 

galaxies in Illustris and then checked at what redshift the average 
ionized fraction in a sphere of radius 150 kpc around them (taking 
their z - 6 positions) reached 50% for the first time. The results 
of this analysis are shown in Figure 6, binned as a function of stel¬ 
lar mass. There is a clear trend for an earlier reionization around 
more massive galaxies. Interestingly, the spread in the reionization 
times slightly increases towards smaller stellar masses, indicating 
that dwarf galaxies are expected to show larger diversity in their 
reionization histories. Also, these low mass galaxies are more sen¬ 
sitive to the adopted parametrization of the escape fraction and tend 
to reionize later in our V5 model. 


3.2 UV background 

In Figure 7, we consider the time evolution of the volume aver¬ 
aged photoionization rate for models calculated with different es¬ 
cape fractions and grid resolutions, and compare to different mod¬ 
els in the literature for the evolution of the cosmic UV background 
(Haardt & Madau 2012; Faucher-Giguere et al. 2009). At high red- 
shift, the background builds up exponentially with redshift, simi¬ 
lar to the growth of the volume weighted ionization fraction, with 
an overall amplitude that varies with the escape fraction model. 
Our variable escape fraction models and our fixed escape fraction 
model bracket the scenario of HM12. Interestingly, our scenario 
V5_1024_M1 follows the model of Faucher-Giguere et al. (2009) 
very closely, except that after reionization is completed, our cal¬ 
culations tend to overshot the model predictions. We note that a 
very similar behaviour is also seen in the recent radiative transfer 
simulations of Pawlik et al. (2015), where this effect is even more 
pronounced. A more realistic variable escape fraction model than 
the rather simple parametrization employed might help resolving 
this issue. 

3.3 Optical depth r 

Starting with the onset of reionization, CMB photons will Thomson 
scatter off the free electrons again. This effect can be quantified by 



Figure 7. Average photoionization rate T as function of redshift for models 
calculated with different escape fractions and grid resolutions. For com¬ 
parison, we also include two widely used theoretical models for the UV 
background evolution: FG09 (Faucher-Giguere et al. 2009, we show the 
updated version from Dec. 2011) and HM12 (Haardt & Madau 2012). The 
V5_1024_M1 model is in particularly good agreement with FG09, at least 
for z < 10. 



0 5 10 15 20 

z 

Figure 8. Cumulative optical depth for Thompson scattering on free elec¬ 
trons, integrated out to the redshift specified on the horizontal axis. Solid 
lines include electrons from doubly ionized helium (assuming that they con¬ 
tribute for z < 3), while dotted lines assume hydrogen and one helium elec¬ 
tron only. The horizontal line with +lcr uncertainty region (shaded) marks 
the newest 2015 constraints r = 0.066 + 0.016 by the Planck Collaboration 
et al. (2015). Out fiducial model V1_1024_M1 is in very good agreement 
with optical depth inferred from these precision measurements of the CMB. 
Our other models lie slightly lower, however their value is still consistent. 
Interestingly, the previous determination by Planck based on their first 2013 
data analysis had given a considerably higher value for r. The tension with 
this result is now resolved. 
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measuring the cumulative optical depth r seen by CMB photons 
along their path towards us. This optical depth is given by 


T = CO-th 


f 


ne(z)—dz, 

dz 


(15) 


where (Xth is the Thomson cross section, and the number density 
of free electrons. In this work, we only consider hydrogen reion¬ 
ization for simplicity. Given that the first ionization potential of 
helium is very close to the ionization potential of hydrogen, we as¬ 
sume that Hell is created in proportion to HIT The free electron 
density is then given by 


He 


1 + 


1-X 

4X 


(1 - vm) wr ~ 1.079 (1 - vm) wr, 


(16) 


where X = 0.76 is the hydrogen mass fraction. At late times, 
helium reionization will eventually be completed, increasing the 
optical depth slightly compared to the above estimate. Adopting 
the common assumption that helium becomes doubly ionized at 
z = 3 (e.g. Iliev et al. 2005), the free electron density increases by 
Ane = (1 - X)/(4X) (1 - xm) nu ^ 0.079 (1 - xu) wr compared to 
the estimate above, and hence the optical is enlarged by 

Ar^ccTth r AfZg(z)^dz = 0.0011, (17) 

Js dz 

which is a very small correction given the other uncertainties. 

The optical depth r for a specific reionization history can be 
converted into an effective reionization redshift Zreion assuming a 
fiducial scenario in which the reionization transition is instanta¬ 
neous at this epoch. The latest WMAP9 results find a best-fit value 
of T = 0.088 + 0.013, corresponding to Zreion = 10.5 + 1.1 (Hinshaw 
et al. 2013), quite a bit lower than the value of r = 0.17 + 0.08 
WMAPl had initially estimated. The results of the PLANCK mis¬ 
sion in its 2013 data release (Planck Collaboration et al. 2014) 
favour a very similar, slightly larger value for the optical depth, 
T = 0.089 + 0.032, corresponding to an even earlier reionization 
redshift of Zreion = 10.8. However, the latest PLANCK data release 
of 2015 (Planck Collaboration et al. 2015) prefers a much lower 
optical depth of r = 0.066 + 0.016 and a corresponding redshift 
of reionization of Zreion = 8 . 8 ^} 4 . A similar low optical depth of 
T = 0.063 + 0.013 has been found in Finkelstein et al. (2014) based 
on a UV luminosity function derived from Hubble Ultra Deep Field 
and Hubble Frontier Field data. 

In Figure 8 , we show the optical depth of our reionization 
simulations as a function of the integration redshift z for three 
of our models. The most recent 2015 constraint from PLANCK is 
shown as a horizontal line, together with the +lcr uncertainty re¬ 
gion (shaded). Our models with a variable escape fraction are com¬ 
fortably compatible within the error bars with the 2015 Planck re¬ 
sults. Our fiducial model V1_1024_M1 predicts an optical depth of 
T = 0.065 which is in very good agreement with the most recent 
Planck 2015 data. However all of our other models prefer the low 
side of the range determined by Planck. Still, it is very promising 
that the former tension between galaxy formation simulations and 
optical depth inferred from CMB measurements seems nearly re¬ 
solved with the 2015 Planck data. A similar finding has been re¬ 
ported in Robertson et al. (2015) based on Hubble observations 
of distant galaxies. Allowing for additional high redshift star for¬ 
mation could easily close the small remaining gap if needed, but 
we note that constraints from galaxy formation disfavour this solu¬ 
tion. For example, Illustris already tends to overshoot estimates for 
the stellar mass function of small galaxies, at late and early times 
alike (Vogelsberger et al. 2014a; Genel et al. 2014). Resolving these 



Figure 9. Distribution of the characteristics radius of ionized regions at 
four different redshifts, in our fiducial model V1_1024_ML The integral 
over the distributions is normalized to unity for each of the measurements, 
and a horizontal line corresponds to equal volume fraction per logarithmic 
size interval. Initially, small bubbles dominate but over time the distribution 
shifts to ever larger bubbles until only one large region dominates. 



Figure 10. The same information as in Fig. 9, except with a different nor¬ 
malization. Here the integral is normalized by the constant comoving box 
volume so that the area under the distributions gives the total ionized vol¬ 
ume fraction. This is informative because the relative constancy of the size 
distributions for small bubble sizes suggests that the growing bubbles are re¬ 
plenished by new small bubbles just at the right rate to achieve this balance. 
Interestingly, the volume occupied by small bubbles stays hence roughly 
the same for the whole duration of reionization, even slightly beyond bub¬ 
ble percolation where a dominating single large ionized region forms. 


problems seems to call for reduced high redshift star formation and 
not the opposite, highlighting the difficulty to reconcile high optical 
depths values from CMB experiments with detailed galaxy forma¬ 
tion models. 


3.4 Bubble size statistics 

Mapping the epoch of reionization more directly than possible thus 
far, for example through 21 cm imaging, is an exciting observational 
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prospect. Once this becomes possible with future radio telescopes 
such as the SKA, quantitative measures of the geometry of the ion¬ 
ized regions, such as their topology, promise to be a powerful probe 
of theoretical models for reionization. Radiation transfer models 
like those calculated herein are the method of choice to make the 
required detailed predictions about how these morphological mea¬ 
sures evolve in time. To illustrate this, we here compute a few basic 
statistical measures that quantify the number and size of ionized 
regions as a function of time, which may also serve as a useful 
comparison against other theoretical reionization models. 

We start by tagging cells as ionized if their ionization fraction 
exceeds 1 - xm > 0.5. Based on the resulting grid of binary values, 
ionized regions are then identified using a friends-of-friends algo¬ 
rithm that links adjacent ionized cells (Iliev et al. 2006; Chardin 
et al. 2012). We let cells belong to the same group if they share 
at least one corner, or in other words, if any of the 26 neighbours 
of an ionized cell is also ionized, both cells are put into the same 
group. For each ionized island identified in this way, we compute 
an effective radius as r = [3/ (4;rV)]^^^, where V is the cumulative 
volume of the cells making up the group. Finally, we consider the 
distribution function fy of the ionized volume fraction contained in 
regions of a given bubble radius. 

In Figure 9, we show the resulting distribution function when 
the convention of Zahn et al. (2007) is followed and the distribution 
is normalized such that 

J' /vrdlogr = 1, (18) 

i.e. we only consider the ionized volume of the box. The results 
show that at early times, when only a small fraction of the volume 
is ionized, the ionized volume is comprised of disjoint regions of 
characteristic size r ^ 2 cMpc. While reionization progresses, ever 
larger bubbles appear, with a flat distribution as a function of size, 
i.e. roughly the same amount of ionized volume is contained per 
logarithmic interval in bubble size, up to bubble sizes of the order 
of r ^ 20 cMpc. Eventually, the bubbles start to percolate and one 
dominating region containing a substantial fraction of the simula¬ 
tion volume is formed (z - 8.2). There is then still a population 
of smaller ionized regions left, with a constant volume fraction per 
unit log r over a dynamic range of about ~ 5 in size. 

An alternative normalization for the size distribution is used 
in Gnedin & Kaurov (2014), who consider the whole box thus that 
integrating over fy gives the ionized volume fraction: 


/ 


/;rdlogr= 1 -xy. 


(19) 


It is instructive to plot the corresponding distributions also with this 
normalization, which is shown in Figure 10. Now the area under the 
distribution function grows with redshift, reflecting the increase of 
the ionized volume fraction. Interestingly, we see in this represen¬ 
tation that the bubble size distribution is fairly constant with time, 
especially for the small bubble sizes. Even though these bubbles 
grow individually in size with time, the fact that the abundance of 
bubbles of a given size stays approximately constant in time (once 
the first bubbles of this size have formed), suggests that small bub¬ 
bles are reformed just at the right rate to compensate for the loss of 
bubbles of a given size due to the growth or coalescence of bubbles. 


3.5 Distribution function of neutral volume fraction 

In Eigure 11, we show the distribution function f{xy) - dV/dvv 
of the volume fraction that is found at a given neutral fraction. 




Figure 11. Degree of ionization level PDFs at different redshifts. The top 
panel shows the results for the V1_1024_M1 model, whereas the bottom 
panel is for the V5_1024_M1 model. Compare to Fig. 7, where the VI run 
shows a strong upturn in the UV flux, which is here expressed as a shift of 
ionized cells down to nni = 10“^ at around z = 5.3, which does not happen 
in the V5 run. 


We compare our two different variable reionization scenarios, in 
the form of V1_1024_M1 (top panel) and V5_1024_M1 (bottom 
panel), and give results for different redshifts in each case. Note that 
we plot vy f{xy) on the vertical axis versus the log of vy, i.e. the 
area under the curves is proportional to the volume fraction at the 
corresponding range of neutral fractions. 

Progress in reionization is associated with a large increase in 
the volume fraction found at low neutral fractions, as is of course 
expected. Interestingly, the differential distribution of volume at a 
given neutral fraction is however fairly broad while reionization 
is not completed, with a peak at a characteristic neutral fraction 
that shifts progressively to lower values. For example, at redshifts 
z ^ 11, most of the volume is either still neutral or at a neutral 
fraction of vy ~ 10"^. At around redshift z ~ 7, the characteris¬ 
tic neutral fraction where most of the volume is found has dropped 
to vy ~ 10""^. Finally, post reionization, the two models start to 
differ more prominently. Here V1_1024_M1 shows a low neutral 
fraction of Xy ~ 10“^ or lower for most of its volume, which corre¬ 
sponds also to the strong rise in the predicted UV background seen 
in this model in Figure 7. In contrast, the model V5_1024_M1, 
which shows good agreement with the UV background model of 


MNRAS 000, 1-18 (2015) 




















Hydrogen Reionization in Illustris 13 


Faucher-Giguere et al. (2009) at this epoch, yields a markedly dif¬ 
ferent distribution of the neutral fraction, with most of the volume 
having neutral fractions around xy - 5 X 10"^. 


4 CAVEATS AND DISCUSSION 

The accuracy and reliability of our results are influenced by many 
numerical aspects as well as physical uncertainties. In the follow¬ 
ing we discuss a number of these aspects, focusing primarily on 
those pertaining to the radiative transfer modelling itself. We note 
however that there are in principle additional uncertainties related, 
for example, to the treatment of star formation and the associated 
feedback processes in the underlying Illustris simulation, or to the 
cosmological background model that we use. These are arguably 
subdominant compared to the uncertainties related to the reioniza¬ 
tion calculation itself (such as escape fraction, radiative transfer 
solver, etc.), and in any case are beyond the scope of this paper (a 
discussion of the uncertainties in the galaxy formation model can 
be found in Vogelsberger et al. 2014a). 

4.1 Reionization feedback 

Due to the fact that we simulate reionization only in post¬ 
processing, any back reaction onto the gas due to photoionization 
heating and potentially radiation pressure is not taken into account 
self-consistently. The Illustris simulation assumes a uniform global 
UV background, hence the average back reaction on star formation 
due to photo-ionization is approximately accounted for, but any lo¬ 
cal modulation of the corresponding effects is of course ignored. 
This limitation could only be overcome by dynamically coupling 
the radiative transfer solver to a hydrodynamical code and doing 
full radiation-hydrodynamics simulations of galaxy formation. Re¬ 
cently, impressive progress has been made in this direction (Gnedin 
2014; Gnedin & Kaurov 2014; Pawlik et al. 2015), but the achieved 
cosmological volumes are still severely limited due to the demand¬ 
ing computational cost of radiative transfer, and in general, these 
calculations have not been evolved to redshift z = 0, and thus it 
is unclear whether they are successfully reproducing the observed 
galaxy population. 

Besides photo-heating, the ionizing radiation of young stars 
could also exert significant feedback effects through radiation pres¬ 
sure, particularly in dusty gas where infrared radiation may be 
trapped (Murray et al. 2010; Hopkins et al. 2012; Agertz et al. 
2013). However, the effectiveness of this mechanism is debated, 
with a number of recent studies arguing that photo-heating is likely 
the dominant feedback channel on the scale of galaxies, with radi¬ 
ation pressure being comparatively unimportant (Sales et al. 2014; 
Rosdahl et al. 2015). We thus consider the omission of radiation 
pressure effects in our reionization calculations to be comparatively 
unimportant. 

4.2 Moment-based versus cone-based RT method 

Our radiative transfer implementation supports two different trans¬ 
port methods, allowing us to compare them against each other with 
no changes in any other aspect of the modelling. In the case of the 
cone-based method, we have to store and process 48 radiation in¬ 
tensity fields, one for each advection direction. On the other hand, 
the moment-based Ml scheme only requires 4 fields, one for the 
photon number density and 3 for the flux vector. This difference 
makes the cone-based advection scheme much more expensive in 


Cone-based Ml closure 
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Figure 12. Visual comparison of our two different radiation advection meth¬ 
ods at two instances in time. All panels show a thin slice through the box 
with a side length of 21.3 cMpc. The top row compares the neutral hydrogen 
fraction at z = 8.5 for the cone-based method (left) and the moment-based 
method with Ml closure (right). Similarly, the bottom row shows this com¬ 
parison for redshift z = 8. 


terms of computational cost as well as in terms of (GPU) memory 
requirements. 

If we compare the ionization histories predicted by these dif¬ 
ferent radiative transfer methods in terms of the ionized volume 
fraction, no appreciable differences are detected. In fact, the agree¬ 
ment of the evolution of the ionized volume fractions is so good 
that we refrain from showing the corresponding comparison in a 
dedicated plot. But this consistency is perhaps not too surprising. 
The ionization history mainly depends on the source population 
that injects ionizing photons, as well as on the density evolution 
of the gas, as the latter sensitively determines the recombination 
rate. Given that the photon injection rates and the density structure 
are exactly equal in our comparison, and given the fact that both 
radiative transfer algorithms are manifestly photon conserving, any 
difference between the cone-based and the Ml-closure methods can 
only be induced by differences in the photon transport directions. 
These are apparently subtle enough that they do not matter much 
for global statistics of the reionization transition. 

However, despite this good agreement in global averages, the 
two methods still show differences in the detailed morphology of 
the ionized bubbles when examined in detail. In Figure 12, we com¬ 
pare the morphology of the ionized regions around a typical galaxy 
at redshifts z = 8 and z = 8.5. While there is clearly a great deal 
of similarity, the detailed locations of the ionization fronts differ 
substantially, highlighting that the radiative transfer solvers do not 
behave identically after all. This is also borne out by a higher or¬ 
der quantitative comparison of the neutral hydrogen fraction fields 
predicted by the two methods. This can for example be done by 
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Figure 13. Impact of the reduced speed of light approximation on the evolu¬ 
tion of the neutral and ionized volume fractions. Reducing the speed of light 
by a factor of 10 (blue line) or 100 (red line) relative to our default calcula¬ 
tion with a physical speed of light (green line) causes a later reionization of 
the universe. 

computing the mass-weighted standard deviation of the difference 
between the neutral hydrogen densities obtained by our two radia¬ 
tive transfer schemes: 

V = var (p (xmi - J^cone)) / <P> • (20) 

For our Vl_0512 models, we find that this quantity rises with de¬ 
creasing redshift until a maximum of v = 0.2 is reached at z = 7.5. 
Afterwards, the full volume is quickly reionized and the variance 
of the difference field rapidly declines again. We note that the cone- 
based method should be the more accurate approach in this com¬ 
parison, as it can avoid certain inaccuracies of the Ml approach, 
in particular when the ionization bubbles of two or more sources 
overlap. 

4.3 How accurate is the reduced speed of light 
approximation? 

As long as the medium is dense enough, the propagation speed of 
the ionization fronts is determined by the rate at which new ioniz¬ 
ing photons arrive at the edge of neutral gas, and not how fast they 
get there. This motivates the idea of the so-called reduced speed 
of light approximation (Gnedin & Abel 2001; Aubert & Teyssier 
2008), in which the physical value of c is artificially reduced. The 
computational advantage of a reduced speed of light is that a much 
larger Courant time step is allowed in schemes where photon trans¬ 
port is followed with explicit time integration. Rosdahl et al. (2013) 
report that the reduced speed of light approximation describes the 
solution of Stromgren sphere well after an effective crossing time 
Across = ^s/c, where rs is the radius of the corresponding Strom¬ 
gren sphere and c the (reduced) speed of light. Before Across = ^s/c, 
however, the numerical solution necessarily always falls behind the 
correct one. Considering the relevant time scales that have to be 
resolved in cosmic reionization, this yields a criterium for the max¬ 
imum allowed reduction of the speed of light. The conclusion of 
Rosdahl et al. (2013) is that there is not much room for applying 
the reduced speed of light approximation if accurate reionization 
simulations of the IGM are desired. 

It is interesting to use our independent radiative transfer code 
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Figure 14. Clumping factor of the gas computed in different ways. The up¬ 
per solid line shows the intrinsic clumping factor of all non star-forming 
gas (with density below the star formation threshold), based on the Voronoi 
tessellation of the Illustris simulation volume. The dotted line gives the in¬ 
trinsic clumping factor of all the gas in the simulation. For comparison, we 
also show the clumping factor of all the non star-forming gas when this is 
obtained for different mesh resolutions after mapping the simulation volume 
to a grid with fixed spatial resolution. 

to check this assessment. In Figure 13, we show the impact of the 
reduced speed of light on the obtained reionization history, based 
on a reduction of the physical speed of light by a factor of 10 or 
100, respectively. Consistently with the findings of Rosdahl et al. 
(2013), the reduction of the speed of light leads to a significant 
delay in the resulting epoch of reionization, amounting to Az ~ 
1-1.5 for the factor of 10 reduction, and much larger for a factor of 
100. The size of this error unfortunately implies that this numerical 
trick can induce unacceptably large distortions in the reionization 
predictions, hence we have refrained from using it throughout the 
study. 

4.4 Clumping factors 

The degree of gas clumping is a critical factor in models of reion¬ 
ization, as it directly determines the recombination rate and hence 
the amount of photons required to reionize the universe and to keep 
it ionized. Full hydrodynamical simulations of galaxy formation are 
a particularly powerful tool to realistically predict the non-linear 
density structure of the gas, and hence to quantify the gas clump¬ 
ing. In Figure 14, we show the full clumping factor C of all the gas, 
defined in the standard way as 


where the averages are volume averages for the full simulation 
box. The black lines represent the clumping factor obtained directly 
from the Voronoi cells of the underlying Illustris simulation, which 
is hence accounting for all structure resolved by the more than 6 bil¬ 
lion cells of the simulation. We give results for the complete density 
field (dotted line) as well as for cells constrained to not lie on the 
effective equation of state of star forming gas (solid black line). The 
former result includes the collapsed gas that corresponds to the ISM 
and is star-forming, whereas in the latter this phase is excluded. The 
clumping factor is obviously much higher when this star-forming 
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Figure 15. Gas clumping factors measured for gas below a certain characteristic gas overdensity. The left-hand panel shows C 20 , with the solid black line giving 
the result for the actual Voronoi tessellation used in the Illustris simulation for the hydrodynamical calculations. The other solid lines give the corresponding 
result for the binned density field when grids with resolution from 256^ to 1024^ are used, as labeled. The right-hand panel gives the same results for Cioo, 
where instead a density threshold of 100 times the mean baryonic density is adopted. For reference, we also show fitting models by Finlator et al. (2012, FI2) 
and (Pawlik et al. 2009, P09), which give the clumping factor of ionized gas or gas below a overdensity threshold of 100 (their z ~ 7.5 reionization case), 
respectively. 


gas is included, but as the sub-grid model used by Illustris glosses 
over the true multi-phase nature of the ISM, the resulting clump¬ 
ing factor for all the gas is still an underestimate. However, we are 
here really only interested in the non-starforming gas, because the 
recombinations and absorptions inside the star forming regions are 
collectively accounted for by the escape fraction, which in part may 
be viewed as parameterizing our ignorance of the detailed gas struc¬ 
ture on ISM scales. 

The solid red, green and blue lines in Fig. 14 show the clump¬ 
ing factor of all gas after binning it on radiation transfer grids with 
different resolution. Clearly, even for the 1024^ grid we loose al¬ 
most a factor of two in the total clumping due to the smoothing of 
this grid. However, as cosmic reionization is a volume filling pro¬ 
cess and the densest gas occupies only a tiny fraction of the value, it 
makes more sense to refer the escape fraction to a somewhat lower 
density threshold than the star-formation threshold. We are hence 
really interested in the clumping factor of gas up to some limited 
overdensity value, for example up to 20 or 100 times the mean bary¬ 
onic density. Since both of these fiducial values have been used in 
the literature, we show in Figure 15 our results for C 20 and Cioo (in 
the left- and right-hand panels, respectively), where only gas cells 
which have a density of at most 20 x pb or 100 x pb have been in¬ 
cluded, respectively, with pb = i^bPcrit denoting the mean baryon 
density. As before, we show results for the underlying Voronoi tes¬ 
sellation as well as reionization grids between 256^ and 1024^ res¬ 
olutions. Note that in these plots a linear scale for the clumping 
factor has been used. 

Our high resolution radiative transfer grid underpredicts the 
C 20 clumping a bit due to the smoothing effects of the binning, but 
the effect is minor. The situation is a bit worse if one also wants to 
get the correct clumping factor for gas up to an overdensity of 100. 


Here some of this additional clumping is not resolved by the 1024^ 
grid, as the results in the right-hand panel of Fig. 15 show. Given the 
trend with increasing resolution, using a 2048^ grid instead (which 
we unfortunately cannot carry out due to memory constraints on 
the GPU system we have presently access to) should however be 
able to fully recover the Cioo clumping of this gas. As we discuss 
in more detail below, this resolution problem for the Cioo quantity 
affects cosmic reionization however only mildly and is hence com¬ 
paratively benign. 

It is interesting to compare our clumping factors with those 
inferred from other works. Finlator et al. (2012) have pointed out 
that different definitions of the clumping factor can produce sub¬ 
stantial differences in the results. It is thus important to base any 
such comparison on the same definition, which sometimes corre¬ 
sponds to considering the clumping of all the gas below a certain 
density threshold (to separate collapsed and diffuse gas), or to re¬ 
stricting the evaluation of the clumping factor to ionized gas. Note 
that the latter depends both on the detailed reionization model and 
the gas distribution. To get a sense of how well the gas distributions 
compare, it is thus arguably best to compare the total gas clumping 
factor. For the Illustris simulation at z = 8, we measure for the 
clumping of the non-starforming gas 10.2, slightly higher than the 
value reported by Finlator et al. (2012). However, our value is sig¬ 
nificantly higher than the total gas clumping factor of C ~ 3 at 
z = 8 reported by Jeeson-Daniel et al. (2014), which makes it con¬ 
siderably easier to achieve reionization in their model. When the 
baryon density is restricted to lie below an overdensity of 100, we 
find a clumping factor of about 4 at z = 8, somewhat larger than 
what was found in Pawlik et al. (2009) for their models reionizing 
at or before z = 9, but a bit lower than their model reionizing at 
z = 7.5. For a yet lower overdensity threshold of 20, Wise et al. 
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(2014) report a value of 6.5, which is above our measurement of 
^ 2.4 for C 20 at redshift z = 8. This likely reflects the higher mass 
resolution of their simulation, which has a boxsize of just 1 Mpc, 
but it could also be affected by the different feedback models in the 
two simulations. 


4.5 Spatial resolution and convergence 

Because the recombination rate depends nonlinearly on the density 
in a cell, the smoothing of density fluctuations (for example as a 
result of binning) causes an underestimate of recombination events 
and hence biases reionization towards higher redshift. Our results 
for the clumping factor indicate that our radiative transfer calcu¬ 
lations clearly suffer from this effect to some degree. However, it 
is not obvious whether the size of the bias is quantitatively signif¬ 
icant in the end, because the clumping of the volume-filling gas 
(which has comparatively low overdensity) is captured well by our 
high-resolution grid. If reionization would mostly occur ‘outside- 
in’, with low density regions ionized early, one may hope that this is 
already sufficient for allowing converged predictions of the reion¬ 
ization redshift even if density peaks are washed out. However, 
given that our results have confirmed that dense regions tend to 
be ionized earlier, this may largely be wishful thinking. 

Indeed, this is borne out by our convergence tests for the reion¬ 
ization history of models V1 and V5 shown in Figure 16. Evidently, 
as the resolution of the grid for the radiative transfer simulation is 
increased, reionization happens progressively somewhat later, as a 
result of the smaller degree of suppression of the true underlying 
dumpiness of the gas. This prevents us from achieving a formal 
numerical convergence for our reionization histories. We note how¬ 
ever that full radiation hydrodynamics simulations will counter this 
drift by physically inducing a reduction of the dumpiness of the 
gas (Pawlik et al. 2009), due to the photo-heating and the resulting 
pressure smoothing. This effect is not included in our simulations 
prior to redshift z = 6, when reionization happens due to the exter¬ 
nally imposed UV background. We thus expect that our 1024^ and 
512^ grids may well bracket the true behavior of the dumping in a 
self-consistent simulation with full radiation hydrodynamics. This 
then also means that a 2048^ calculation without taking this effect 
into account may well produce a less accurate result than the 1024^ 
grid we used. 


5 CONCLUSIONS 

In this work, we considered only hydrogen reionization and focused 
on ordinary high-redshift star formation as the primary source of 
ionizing photons. Other populations may in principal contribute to 
reionization, in particular primordial population-III stars, AGNs, 
or exotic sources such as annihilating dark matter. While Pop-Ill 
stars may be important for the onset of reionization, most estimates 
for their relative contribution to global star formation suggest that 
neglecting them for reionization is justified. Still, accounting for 
them in future models would be dearly desirable, if only for com¬ 
pleteness. Our neglect of AGN radiation for hydrogen reionization 
is however quite well justified because their ionizing luminosity is 
overwhelmed by star formation at high redshift. AGNs become im¬ 
portant at intermediate redshifts, however, where they likely play 
an important role in Hell reionization. 

Even with these simplifying assumptions, calculating the ion¬ 
izing flux that becomes available for reionization is affected by a 



z 


Figure 16. Convergence study for the neutral and ionized volume fractions 
for our reionization simulations based on the Ml method. The green and 
red lines are for the two different variable escape fraction models we consid¬ 
ered. There are significant residual trends with resolution due to the smooth¬ 
ing effects a coarse grid has on the clumpy gas distribution. As a result, the 
highest resolution simulation tends to reionize slightly later than predicted 
by calculations at lower resolution. 


substantial number of uncertainties. This includes the stellar popu¬ 
lation synthesis model we have used, and in particular, the adopted 
stellar initial mass function, where we used a Chabrier IMF and as¬ 
sumed that there are no significant IMF variations as a function of 
environment. Another major uncertainty lies in the escape fraction, 
which is physically uncertain and is in large part a phenomeno¬ 
logical parameter in our models, absorbing uncertainties due to the 
treatment of dense gas and the limited spatial resolution. Finally, 
there are also numerical limitations related to the radiative trans¬ 
fer solver, the finite angular and spatial resolution employed in the 
radiative transfer, and the lack of self-consistently accounting for 
local feedback effects by the radiation held. 

Fortunately some of these uncertainties can be greatly reduced 
by matching key observables such as the total optical depth for elec¬ 
tron scattering or the amplitude of the metagalactic ionizing UV 
background after completion of reionization (e.g. Faucher-Giguere 
et al. 2008a,b). Our radiative transfer models for Illustris then basi¬ 
cally test whether cosmic reionization does occur for reasonable as¬ 
sumptions about the escape fraction, based on a galaxy population 
that yields a successful description of a slew of other observational 
data both at high and low redshift. To the extent this works, it pro¬ 
vides reassurance for the cosmological consistency of our galaxy 
formation and reionization models, and it shows that they are physi¬ 
cally meaningful. Importantly, they can hence be used to learn more 
about the reionization transition itself. Our main findings can be 
summarized as follows: 

(i) The star formation history predicted by the Illustris simula¬ 
tion combined with its high-resolution gas density held allows cos¬ 
mic reionization with an optical depth of r = 0.065, consistent 
with the latest Planck 2015 results as well as constraints from high- 
redshift quasars. This relies on ordinary stellar populations only, 
but requires optimistic assumptions for a high escape fraction at 
high redshift. 

(ii) Previous tensions between the high optical depth favoured 
by CMB results and the low level of high redshift star formation 
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required by successful galaxy formation models are thus essentially 
resolved with Planck 2015. 

(iii) Using a suitable variable escape fraction model, we can ap¬ 
proximately reproduce the expected UV background after reion¬ 
ization is completed, but most of our models tend to then over¬ 
shoot the ionization rate and yield a slightly lower neutral fraction 
than inferred from quasar absorption lines. A fine-tuned model may 
plausibly yield an improved match. 

(iv) Reionization proceeds inside out in our model, with over- 
dense regions being ionized earlier on average than lower density 
regions. 

(v) The size distribution of ionized regions shows a remarkable 
constancy with time for small bubble sizes, suggesting that during 
reionization new small bubbles are formed roughly at the rate at 
which they are removed by size growth or coalescence with other 
regions. The characteristic size of bubbles grows with time, but 
there is a fairly flat distribution of sizes between r ~ 2 cMpc and 
r ~ 20 cMpc with about equal volume fraction per logarithmic size 
interval just before reionization is completed. 

(vi) The duration of the reionization transition varies with the 
escape fraction model, and we find transition times for a drop of 
the neutral fraction from 80 to 20% between 190 and 340 Myr, 
depending on the model. 

(vii) The distribution of volume with respect to neutral fraction 
is quite broad but shows a peak that progresses to ever smaller neu¬ 
tral fraction with time. For models that successfully match the UV 
background constraints after reionization, the characteristic neutral 
fraction is 5 X 10“^, with the lowest amount of volume found at 
neutral fractions of 2 X 10“^. 

In future work based on our methodology it would be particu¬ 
larly interesting to study Hell reionization. Due to the two ioniza¬ 
tion levels of helium with different ionization thresholds, the com¬ 
putational cost rises by at least a factor of about three, as more 
spectral bins have to be tracked instead of just one for hydrogen. 
Additionally, a much longer physical time down to z ~ 3 has to 
be followed. This would still be possible in postprocessing with 
a future Illustris type simulation with a larger box size, contain¬ 
ing an evolving AGN population including contributions from the 
brightest objects. Simulating this much physical time is extremely 
challenging for direct radiation-hydrodynamics simulations, much 
more so than hydrogen reionization simulations that can be stopped 
at z - 6. This makes postprocessing approaches the only practi¬ 
cal radiative transfer method to study Hell reionization in the near 
term. 
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